## describe homophily wrt different variables
describe_sorting<-function(A_in) function(measure)
{
  own_vec<-measure
  n_friends_1<-rowSums(A_in[,,1])
  total_friend_vec_1<- A_in[,,1] %*% own_vec
  mean_friend_vec_1<- total_friend_vec_1/ n_friends_1
  n_friends_2<-rowSums(A_in[,,2])
  total_friend_vec_2<- A_in[,,2] %*% own_vec
  mean_friend_vec_2<- total_friend_vec_2/ n_friends_2
  
  mean_friend_vec<-c(mean_friend_vec_1, mean_friend_vec_2)
  
  cor(c(own_vec,own_vec), mean_friend_vec, use="complete.obs")
}

## make graphs
A_1_graph<-graph.adjacency(A_array_months[,,1])
A_2_graph<-graph.adjacency(A_array_months[,,2])

if(simulate_A_array==F | sim_no_friends==F){
  A_summary_table<-rbind(frac_0_friends=c(obs=mean(apply(A_array_months_data, c(2,3), sum)==0), sim=mean(apply(A_array_months, c(2,3), sum)==0)),            mean_friends=c(obs=mean(apply(A_array_months_data, c(2,3), sum)), sim=mean(apply(A_array_months, c(2,3), sum))),                     median_friends=c(obs=median(apply(A_array_months_data, c(2,3), sum)), sim=median(apply(A_array_months, c(2,3), sum))),               mean_triangles=c(mean(c(adjacent.triangles(graph.adjacency(A_array_months_data[,,1])), adjacent.triangles(graph.adjacency(A_array_months_data[,,2])))), mean(c(adjacent.triangles(graph.adjacency(A_array_months[,,1])), adjacent.triangles(graph.adjacency(A_array_months[,,2]))))), 
                         median_triangles=c(median(c(adjacent.triangles(graph.adjacency(A_array_months_data[,,1])), adjacent.triangles(graph.adjacency(A_array_months_data[,,2])))), median(c(adjacent.triangles(graph.adjacency(A_array_months[,,1])), adjacent.triangles(graph.adjacency(A_array_months[,,2]))))), 
                         mean_degree=c(mean(c(degree(graph.adjacency(A_array_months_data[,,1])), degree(graph.adjacency(A_array_months_data[,,2])))), mean(c(degree(graph.adjacency(A_array_months[,,1])), degree(graph.adjacency(A_array_months[,,2]))))), 
                         median_degree=c(median(c(degree(graph.adjacency(A_array_months_data[,,1])), degree(graph.adjacency(A_array_months_data[,,2])))), median(c(degree(graph.adjacency(A_array_months[,,1])), degree(graph.adjacency(A_array_months[,,2]))))), 
                         cor_black_i_not_i= c(describe_sorting(A_array_months_data)(student_data$black), describe_sorting(A_array_months)(student_data$black)),             cor_male_i_not_i= c(describe_sorting(A_array_months_data)(student_data$male), describe_sorting(A_array_months)(student_data$male)),                      cor_hsgpa_i_not_i= c(describe_sorting(A_array_months_data)(student_data$hsgpa), describe_sorting(A_array_months)(student_data$hsgpa)),             cor_combact_i_not_i= c(describe_sorting(A_array_months_data)(student_data$combact), describe_sorting(A_array_months)(student_data$combact)), 
                         cor_studyhs_i_not_i= c(describe_sorting(A_array_months_data)(student_data$studyhs), describe_sorting(A_array_months)(student_data$studyhs)), 
                         cor_estudy_i_not_i= c(describe_sorting(A_array_months_data)(student_data$estudy), describe_sorting(A_array_months)(student_data$estudy)), 
                         cor_mu_si_not_i= c(describe_sorting(A_array_months_data)(solution$mu_study_model), describe_sorting(A_array_months)(solution$mu_study_model)))
  
  print(A_summary_table)
  print(xtable(A_summary_table, caption="Properties of observed and simulated networks"), file="./output/A_summary_table.tex")
  write.csv(xtable(A_summary_table, caption="Properties of observed and simulated networks"), file="./output/A_summary_table.csv")
}

avg_num_friends<-(apply(A_array_months_data, 1, sum)/2)
closeness_A_1<-closeness(A_1_graph)
closeness_A_2<-closeness(A_2_graph)
avg_closeness<-(closeness_A_1 + closeness_A_2)/2

summary(lm(avg_num_friends ~ male + black + hsgpa + combact + studyhs + estudy, data=student_data, na.action="na.exclude"))
summary(lm(avg_closeness ~ male + black + hsgpa + combact + studyhs + estudy, data=student_data, na.action="na.exclude"))
summary(lm(avg_num_friends ~ solution$mu_study_model + I(solution$mu_study_model^2)))
summary(lm(avg_closeness ~ solution$mu_study_model + I(solution$mu_study_model^2)))
summary(lm(closeness_A_1 ~ solution$mu_study_model + I(solution$mu_study_model^2)))
summary(lm(closeness_A_2 ~ solution$mu_study_model + I(solution$mu_study_model^2)))
